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Abstract 


Decision  rules  for  segmenting  scenes  and  for  detecting  the  presence  of  distin¬ 
guished  objects  in  digital  images  can  be  based  on  classical  principles  of  statistical 
inference  if  appropriate  mathematical  models  are  developed  for  observable  pictures. 
The  main  goal  of  this  research  was  to  devise  and  analyze  alternative  image  mod¬ 
els  for  digitized  FLIR  images.  The  work  has  been  done  in  close  cooperation  with 
the  Advanced  Modeling  Team  of  the  U.S.  Army  Night  Vision  and  Electro-Optics 
Laboratory,  Ft.  Belvoir,  Virginia.  This  report  concentrates  on  hierarchical  Markov 
Random  Field  models  and  their  application  to  restoration  and  segmentation  of  FLIR 
images. 
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1  INTRODUCTION. 

Our  primary  goal  has  been  to  construct  a  mathematical  foundation  for  the  rational 
choice  of  decision  functions  for  image  analysis.  This  has  included  structured  models 
for  the  background  against  which  certain  objects,  such  as  tanks,  trucks,  or  armored 
personnel  carriers,  appear.  The  backgrounds  are  “complex”  in  that  their  composition 
is  highly  variable  and  cannot  be  known  in  advance.  The  objects  are  “simple”  in  that 
they  can  be  characterized  by  a  small  number  of  parameters.  While  the  emphasis  has 
been  placed  on  the  logical  and  mathematical  foundations, considerable  effort  has  been 
given  to  the  construction  of  algorithms.  It  is  important  to  keep  the  algorithmic  issues  in 
mind  so  that  we  arrive  at  decision  procedures  that  work  and  that  can  be  computed  with 
reasonable  resources. 

This  report  focuses  on  a  strategy  for  image  modeling  that  has  been  developed  for 
a  number  of  practical  settings.  Here  we  develop  it  for  the  analysis  of  FLIR  images. 
Indeed,  this  project — while  it  is  immediately  concerned  with  problems  suggested  by  the 
U.S.  Army  Night  Vision  and  Electro-Optics  Laboratory — has  had  a  tremendous  impact 
on  the  development  of  a  general  Bayesian  methodology  for  automatic  analysis  of  digital 
images.  Today  that  methodology  is  successfully  addressing  practical  problems  in  medical 
imaging  (computed  tomography,  ultrasound),  remote  sensing  (interpretation  of  SAR 
images),  automatic  inspection  (analysis  of  textured  optical  images  of  silicon  wafers),  and 
image  understanding  (optical  character  recognition,  boundary  finding,  segmentation). 

In  the  interest  of  presenting  a  self-contained  and  coherent  report  on  mathematical 
models  for  FLIR  images,  we  shall  concentrate  this  paper  on  the  general  Bayesian  model 
and  its  adaptation  to  FLIR  imagery.  Our  interactions  with  the  Advanced  Modeling  Team 
at  NV&EOL  have  had  many  other  facets,  including  frequent  on-site  working  sessions, 
supervision  of  the  development  of  computer  algorithms,  direction  for  the  formation  of  a 
data  base  of  features  of  FLIR  images,  statistical  analyses,  and  assistance  with  providing 
information  on  other  mathematical  modeling  efforts.  These  interactions  are  all  directly 
related  to  the  overall  project  on  image  modeling,  and  are  documented  elsewhere.  In 
particular,  the  internal  working  memoranda  listed  in  Appendix  A  provide  additional 
details  on  both  theoretical  and  practical  aspects  of  the  effort. 

Section  2  of  this  paper  gives  an  overview  and  basic  examples  of  the  Bayesian  modeling 
strategy.  It  covers  the  range  of  issues  from  specification  of  the  probabilistic  framework 
to  the  design  of  computational  algorithms. 

Section  3  describes  the  adaptation  of  the  general  Bayesian  paradigm  to  digitized  FLIR 
images.  Here  we  describe  a  specific  heierarchical  probabilistic  model  which  is  useful  for 
FLIR  image  restoration  and  segmentation. 

Section  4  presents  a  FORTRAN  implementation  of  the  image  restoration  algorithm. 
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Program  listings  are  included. 

Section  5  briefly  describes  the  application  of  the  restoration  algorithm  to  eight  ex¬ 
amples  of  FLIR  images  provided  to  us  by  NV&EOL. 

Finally,  two  appendices  include,  respectively,  (i)  a  list  of  internal  working  papers 
developed  during  the  project  and  previously  shared  with  the  Advanced  Modeling  Team 
at  NV&EOL  and  (ii)  pictures  illustrating  the  examples  cited  in  Section  5. 

We  gratefully  acknowedge  the  contributions  made  to  this  research  effort  by  Frank 
Shields  and  Vince  Mirelli  of  the  Advanced  Modeling  Team  at  NV&EOL.  The  discussions 
of  the  fundamental  mathematical  issues  with  Dr.  Mirelli  have  provided  a  tremendous 
stimulus  for  focusing  our  efforts  on  meaningful  ways  of  bringing  mathematics  to  bear  on 
challenging  practical  problems. 


2  BAYESIAN  PARADIGM 


In  real  scenes,  neighboring  pixels  typically  have  similar  intensities,  boundaries  are  usually 
smooth  and  often  straight,  textures,  although  sometimes  random  locally,  define  spatially 
homogeneous  regions,  and  objects,  such  as  grass,  tree  trunks,  branches  and  leaves,  have 
preferred  relations  and  orientations.  Our  approach  to  picture  processing  is  to  articulate 
such  regularities  mathematically,  and  then  to  exploit  them  in  a  statistical  framework 
to  make  inferences.  The  regularities  are  rarely  deterministic;  instead,  they  describe 
correlations  and  likelihoods.  This  leads  us  to  the  Bayesian  formulation,  in  which  prior 
expectations  are  formally  represented  by  a  probability  distribution.  Thus  we  design 
a  distribution  (a  “prior”)  on  relevant  scene  attributes  to  capture  the  tendencies  and 
constraints  that  characterize  the  scenes  of  interest.  Picture  processing  is  then  guided 
by  this  prior  distribution,  which,  if  properly  conceived,  enormously  limits  the  plausible 
restorations  and  interpretations. 

The  approach  involves  five  steps,  which  we  shall  briefly  review  here  (see  [4]  and 
[9]  for  more  details).  This  will  define  the  general  framework,  and  then,  in  the  following 
sections,  we  will  concentrate  on  the  analysis  of  samples  of  FLIR  images,  as  an  illustrative 
application. 

2.1  Image  Models. 

These  are  probability  distributions  on  relevant  image  attributes.  Both  for  reasons  of 
mathematical  and  computational  convenience,  we  use  Markov  random  fields  (MRF)  as 
prior  probability  distributions.  Let  us  suppose  that  we  index  all  of  the  relevant  attributes 
by  the  index  set  S.  S  is  application  specific.  It  typically  includes  indices  for  each  of  the 
pixels  (about  512x512  in  the  usual  video  digitization)  and  may  have  other  indices  for 
such  attributes  as  boundary  elements,  texture  labels,  object  labels  and  so-on.  Associated 
with  each  “site”  s  €  S  is  a  real-valued  random  variable  X,,  representing  the  state  of  the 
corresponding  attribute.  Thus  X,  may  be  the  measured  intensity  at  pixel  s  (typically, 
X,  e  {0,  ...255}),  or  simply  1  or  0  as  a  boundary  element  at  location  s  is  present  or 
absent. 

The  kind  of  knowledge  we  represent  by  the  prior  distribution  is  usually  “local,”  which 
is  to  say  that  we  articulate  regularities  in  terms  of  small  local  collections  of  variables. 
In  the  end,  this  leads  to  a  distribution  on  X  =  {AT4}jes  with  a  more  or  less  “local 
neighborhood  structure”  (again,  we  refer  to  [4]  and  [9j  for  details).  Specifically,  our  priors 
are  Markov  random  fields:  there  exists  a  (symmetric)  neighborhood  relation  G  =  {G,}5€5, 
wherein  G,  C  S  is  the  set  of  neighbors  of  s,  such  that 

II(X,  =  x,\Xr  =  zr,r  E  S,r  ^  s)  =  n(AT,  =  x,\Xr  =  xr,r  G  G,) 


n(a|6)  is  conditional  probability,  and,  by  convention,  s  $  G„.  G  symmetric  means 
s  €  Gr  O  r  G  Ga.  (Here,  we  assume  that  the  range  of  the  random  vector  X  is  discrete; 
there  are  obvious  modifications  for  the  continuous  or  mixed  case.) 

It  is  well-known,  and  very  convenient,  that  a  distribution  IT  defines  a  MRF  on  S  with 
neighborhood  relation  G  if  and  only  if  it  is  Gibbs  with  respect  to  the  same  graph,  (5,  G). 
The  latter  means  that  II  has  the  representation 


(2.1) 

n(x)  =  -e"^ 
z 

where 

(2.2) 

U(x)  =  £  V,(x) 

cec 


C  is  the  collection  of  all  cliques  in  (5,  G)  (collections  of  sites  such  that  every  two  sites 
are  neighbors),  and  Vc(x)  is  a  function  depending  only  on  {x,},6c.  U  is  known  as  the 
“energy,”  and  has  the  intuitive  property  that  the  low  energy  states  are  the  more  likely 
states  under  IT.  The  normalizing  constant,  z,  is  known  as  the  “partition  function”.  The 
Gibbs  distribution  arises  in  statistical  mechanics  as  the  equilibrium  distribution  of  a 
system  with  energy  function  U. 

As  a  simple  example  (too  simple  to  be  of  much  use  for  real  pictures)  suppose  the 
pixel  intensities  are  known,  a  priori,  to  be  one  of  two  levels,  minus  one  (“black”)  or  plus 
one  (“white”).  Let  5  be  the  N  x  N  square  lattice,  and  let  G  be  the  neighborhood  system 
that  corresponds  to  nearest  horizontal  and  vertical  neighbors: 

o  —  o  - 

I  I 

O  -  o  ... 

I  I 

o  -  0  •  •  • 

For  picture  processing,  think  of  N  as  typically  512.  Suppose  that  the  only  relevant 
regularity  is  that  neighboring  pixels  tend  to  have  the  same  intensities.  An  “energy” 
consistent  with  this  regularity  is  the  “Lsing”  potential: 

U{x)  =  -P^x.it  0  >  0 

(».<) 

where  £(#,«)  means  summation  over  all  neighboring  pairs  s,t  G  5.  The  minimum  of  U  is 
achieved  when  x,  =  xt  Vs,t  6  S.  Under  (2.1),  the  likely  pictures  are  therefore  the  ones 


that  respect  our  prior  expectations;  they  segment  into  regions  of  constant  intensities. 
The  larger  (3,  the  larger  the  typical  region.  Later  we  will  discuss  the  issue  of  estimating 
model  parameters  such  as  (3.  (With  energy  (2.2),  II  in  (2.1)  is  called  the  Ising  model. 
It  models  the  equilibrium  distribution  of  the  spin  states  of  the  atoms  in  a  ferromagnet. 
These  spins  tend  to  “line  up,”  and  hence  the  favored  configurations  contain  connected 
regions  of  constant  spins.) 

One  very  good  reason  for  using  MRF  priors  is  their  Gibbs  representations.  Gibbs 
distributions  are  characterized  by  their  energy  functions,  and  these  are  more  convenient 
and  intuitive  for  modelling  than  working  directly  with  probabilities.  See,  for  example, 
[3],  [4],  [5],  [9],  and  [13]  for  many  more  examples,  and  Section  3  below  for  a  more  complex 
and  useful  MRF  model. 

2.2  Degradation  Model. 

The  image  model  is  a  distribution  II(-)  on  the  vector  of  image  attributes  X  =  {X,},6s. 
By  design,  the  components  of  this  vector  contain  all  of  the  relevant  information  for  the 
image  processing  task  at  hand.  Hence,  the  goal  is  to  estimate  X.  This  estimation  will  be 
based  upon  partial  or  corrupted  observations,  and  based  upon  the  prior  information.  In 
emission  tomography,  X  represents  the  spatial  distribution  of  isotope  in  a  target  region 
of  the  body.  What  is  actually  observed  is  a  collection  of  photon  counts  whose  probability 
law  is  Poisson,  with  a  mean  function  that  is  an  attenuated  radon  transform  of  X.  In  the 
texture  labelling  problem,  X  is  the  pixel  intensity  array  and  a  corresponding  array  of 
texture  labels.  Each  label  gives  the  texture  type  of  the  associated  pixel.  The  observation 
is  only  partial:  we  observe  the  pixels,  which  are  just  the  digitized  picture,  but  not  the 
labels.  The  purpose  is  then  to  estimate  the  labels  from  the  picture.  In  a  generic  model 
for  FLIR  images  described  in  Section  3,  X  is  a  hierarchical  model  built  from  the  pixel 
intensity  array  and  from  a  superimposed  array  of  unobservable  edge  elements.  Again,  the 
observation  is  only  partial:  we  observe  the  pixels,  degraded  as  they  are  by  atmospheric 
effects  and  the  sensor,  but  not  the  edge  elements  that  are  combined  to  form  boundaries 
between  objects  and  background.  A  purpose  of  image  segmentation  is  to  estimate  the 
boundaries  from  the  observed  picture. 

The  observations  are  related  to  the  image  process  (X)  by  a  degradation  model  This 
models  the  relation  between  X  and  the  observation  process,  say  Y  =  For  texture 

analysis,  we  will  define  X  =  {Xp ,XL),  where  Xp  is  the  usual  grey-level  pixel  intensity 
process,  and  XL  is  an  associated  array  of  texture  labels.  The  observed  picture  is  just 
Xp ,  and  hence  Y  =  Xp :  the  degradation  is  a  projection.  More  typically,  the  degradation 
involves  a  random  component,  as  in  the  tomography  setting  where  the  observations  are 
Poisson  variables  whose  means  are  related  to  the  image  process  X.  A  more  simple,  and 
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widely  studied  (if  unrealistic)  example  is  additive  “white”  noise.  Let  X  =  {X,},es  be 
just  the  basic  pixel  process.  In  this  case,  T  =  S,  and  for  each  s£  S  we  observe 


r,  =  +  r). 


where,  for  example,  {r?4}4g5  is  Gaussian  with  independent  components,  having  means  0 
and  variances  o2 . 

Formally,  the  degradation  model  is  a  conditional  probability  distribution,  or  density, 
for  Y  given  X:  n(y|x).  If  the  degradation  is  just  added  “white  noise,”  as  in  the  above 
example,  then 

l£i 

For  labelling  textures,  the  degradation  is  deterministic:  n(y|x)  is  concentrated  on  y  =  xp , 
where  x  =  (xp,xL)  has  both  pixel  and  label  components. 


2.3  Posterior  Distribution. 


This  is  the  conditional  distribution  on  the  image  process  X  given  the  observation  process 
Y .  This  “posterior”  or  “ a  posteriori”  distribution  contains  the  information  relevant  to  the 
image  restoration  or  image  analysis  task.  Given  an  observation  Y  =  y,  and  assuming  the 
image  model  (II(x))  and  degradation  model  (II(y|x)),  the  posterior  distribution  reveals 
the  likely  and  unlikely  states  of  the  “true”  (unobserved)  image  X.  Having  constructed  X 
to  contain  all  relevant  image  attributes,  such  as  locations  of  boundaries,  labels  of  objects 
or  textures,  and  so-on,  the  posterior  distribution  comes  to  play  the  fundamental  role  in 
our  approach  to  image  processing. 

The  posterior  distribution  is  easily  derived  from  “Bayes’  rule” 


n(x|y) 


n(yjx)n(x) 

n(y) 


The  denominator,  n(y),  is  difficult  to  evaluate.  It  derives  from  the  prior  and  degrada¬ 
tion  models  by  integration:  n(y)  =  Jz  n(y|x)n(dx),  but  the  formula  is  computationaly 
intractable.  Happily,  our  analysis  of  the  posterior  distribution  will  require  only  ratios, 
not  absolute  probabilities.  Since  y  is  fixed  by  observation,  is  a  constant  that  can  be 
ignored  (see  paragraph  below  on  “computing”). 

As  an  example  we  consider  the  simple  “Ising  model”  prior,  with  observations  cor¬ 
rupted  by  additive  white  noise.  Then 


n(x) 


1 

-ex 

z 


p{-/?Lx.iJ 

(».<) 
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and 

ill 

The  posterior  distribution  is  then 

n(x|y)  =  —ex p{-(3  ]T  x,x(  -  ^(y,  -  x,)2} 

2P  (,,0  2CT  <65 

We  denote  by  zp  the  normalizing  constant  for  the  posterior  distribution.  Of  course,  it 
depends  upon  y,  but  the  latter  is  fixed.  Notice  that  the  posterior  distribution  is  again 
a  MRF.  In  the  case  of  additive  white  noise,  the  neighborhood  system  of  the  posterior 
distribution  is  that  of  the  prior,  and  hence  local.  For  a  wide  class  of  useful  degradation 
models,  including  combinations  of  blur,  added  or  multiplicative  “colored  noise,”  and  a 
variety  of  nonlinear  transformations,  the  posterior  distribution  is  a  MRF  with  a  more 
or  less  local  graph  structure.  This  is  convenient  for  our  computational  schemes,  as  we 
shall  see  shortly.  We  should  note,  however,  that  exceptions  occur.  In  tomography,  for 
example,  the  posterior  distribution  is  associated  with  a  highly  non-local  graph.  This 
situation  incurs  a  high  computational  cost  (see  [5]  for  more  details). 

2.4  MAP  Estimate. 

In  our  framework,  image  processing  amounts  to  choosing  a  particular  image  x,  given  an 
observation  Y  =  y.  A  sensible,  and  suitably-defined  optimal,  choice  is  the  “maximum  a 
posteriori,”  or  “MAP”  estimate: 

choose  x  to  maximize  fl(xjy) 


K 
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The  MAP  estimate  chooses  the  most  likely  x,  given  the  observation.  In  most  applications, 
our  goal  is  to  identify  the  MAP  estimate,  or  a  suitable  approximation.  However,  in  some 
settings  other  estimators  are  more  appropriate.  We  have  found,  for  example,  that  the 
posterior  mean  (/  xll(dx|y))  is  more  effective  for  tomography,  at  least  in  our  experiments. 
Here,  we  concentrate  on  MAP  estimation. 

In  most  applications  we  can  not  hope  to  identify  the  true  maximum  a  posteriori  image 
vector  x.  To  appreciate  the  computational  difficulty,  consider  again  the  Ising  model  with 
added  white  noise: 


(2.3) 
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Zp  (-.0 


1 


2cr2 


J2(y»  -  x>)2} 


>€S 


-•v .  v  . 


* 

I 


This  is  to  be  maximized  over  all  possible  vectors  x  =  £  {  —  1,1}ISL  With 

S  ~  105,  brute  force  approaches  are  intractable;  instead,  we  will  employ  a  Monte 
Carlo  algorithm  which  gives  adequate  approximations. 

Maximizing  (2.3)  amounts  to  minimizing 

up{x)  =  -0Y,  x>Xt  x*)2 

(>.t) 

which  might  be  thought  of  as  the  “posterior  energy”.  (As  with  zp,  the  fixed  observation  y 
is  suppressed  in  the  notation  Up(x).)  More  generally,  we  write  the  posterior  distribution 
as 

(2.4)  —  exp{-Up{x)} 

zp 

and  characterize  the  MAP  estimator  as  the  solution  to  the  problem 

choose  x  to  minimize  Up(x) 

The  utility  of  this  point  of  view  is  that  it  suggests  a  further  analogy  to  statistical  me¬ 
chanics,  and  a  computation  scheme  for  approximating  the  MAP  estimate,  which  we  shall 
now  describe. 

2.5  Computing. 

Pretend  that  (2.4)  is  the  equilibrium  Gibbs  distribution  of  a  real  system.  Recall  that 
MAP  estimation  amounts  to  finding  a  minimal  energy  state.  For  many  physical  systems 
the  low  energy  states  are  the  most  ordered,  and  these  often  have  desirable  properties.  The 
state  of  silicon  suitable  for  wafer  manufacturing,  for  example,  is  a  low  energy  state.  Phys¬ 
ical  chemists  achieve  low  energy  states  by  heating  and  then  slowly  cooling  a  substance. 
This  procedure  is  called  annealing.  Cerny  ( 1 )  and  Kirkpatrick  [12]  suggest  searching  for 
good  minimizers  of  U (•)  by  simulating  the  dynamics  of  annealing,  with  U  playing  the 
role  of  energy  for  an  (imagined)  physical  system.  In  our  image  processing  experiments, 
we  often  use  “simulated  annealing”  to  find  an  approximation  to  the  MAP  estimator. 

Dynamics  are  simulated  by  producing  a  Markov  chain,  X(l),  X(2), ...  with  transi¬ 
tion  probabilities  chosen  so  that  the  equilibrium  distribution  is  the  posterior  (Gibbs) 
distribution  (2.4).  One  way  to  do  this  is  with  the  “Metropolis  algorithm”  [l4h  More 
convenient  for  image  processing  is  a  variation  we  call  stochastic  relaxation.  The  full  story 
can  be  found  in  4  and  9  .  Briefly,  in  stochastic  relaxation  we  choose  a  sequence  of  sites 
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5(1),  3(2), ...  £  S  such  that  each  site  in  S  is  “visited”  infinitely  often.  If  X(t)  =  x,  say, 
then  Xr(t  +  1)  =  xr  Vr  s(t),r  £  S,  and  X,^)(t  +  1)  is  a  sample  from 

n(X>(t)  =  -|  Xr  =  xr,r?s(t)), 

the  conditional  distribution  on  Xs(t)  given  Xr  =  xr  Vr  ^  s{t).  By  the  Markov  property, 
n(xj(()  =  .|XP  =  xr,r  ^  3(0)  =  n(x.(0  =  -\xr  =  xr,r  e  g^) 


where  {G£},es  is  the  posterior  neighborhood  system,  determined  by  the  posterior  energy 
Up(-).  The  prior  distributions  that  we  have  experimented  with  have  mostly  had  local 
neighborhood  systems,  and  usually  the  posterior  neighborhood  system  is  also  more  or 
less  local  as  well.  This  means  that  |Gyfj|  is  small,  and  this  makes  it  relatively  easy  to 
generate,  Monte  Carlo,  X(t  +  1)  from  X(t).  In  fact,  if  ft  is  the  range  of  X,(t),  then 


(2.5) 

where 


n(*.( 


_  .  I  V  —  ^  c  \  _  n(q,  ,(t)X) 
t)  CL  Xr  Xr,  £  Gj(t))  -  nfA 


n(d,»(<)x) 


(ai  i{t)X)r 


a  if  r  =  s(t) 

Xr  if  r  /  s(t) 


Notice  that  (fortunately!)  there  is  no  need  to  compute  the  posterior  partition  function 
zp.  Also,  the  expression  on  the  right  hand  side  of  (2.5)  involves  only  those  potential 
terms  associated  with  cliques  containing  s(t),  since  all  other  terms  are  the  same  in  the 
numerator  and  the  denominator. 

To  simulate  annealing,  we  introduce  an  artificial  “temperature”  into  the  posterior 
distribution: 


nr(x)  = 


ezp{ 


-I'.tz) 


} 


ZP{T) 


As  T  — ►  0,  nT(-)  concentrates  on  low  energy  states  of  i'p.  To  actually  find  these  states, 
we  run  the  stochastic  relaxation  algorithm  while  slowly  lowering  the  temperature.  Thus 
T  =  T(t),  and  T(t)  „  0.  rir^)(-)  replaces  fl(-)  in  computing  the  transition  X(t)  -» 
X[t  +  1).  In  ;4j  we  showed  that,  under  suitable  hypotheses  on  the  sequence  of  site  visits. 

5(l)>  s(2)5 


If  T(t)  >  TTj’T(t)  ♦  0-  then  for  all  c  sufficiently  large  X(t) 

converges  weakly  to  the  distribution  concentrating  uniformly  on 
{x  :  U(x)  =  minvU(y)}. 


More  recently,  our  theorem  has  been  improved  upon  by  many  authors.  In  particular, 
the  smallest  constant  c  which  guarantees  convergence  of  the  annealing  algorithm  to  a 


global  minimum  can  be  specified  in  terms  of  the  energy  function  l'P  (see  6  and  10  ). 
Also,  see  Gidas  7'  for  some  ideas  about  faster  annealing  via  "renormalization  group" 
methods. 

In  the  experiments  with  FLIR  images  to  be  described  here,  MAP  estimates  are  ap¬ 
proximated  by  using  the  annealing  algorithm.  This  involves  Monte  Carlo  computer- 
generation  of  the  sequence  A'(I),  A'(2), ....  terminating  when  the  state  ceases  to  change 
substantially. 


3  A  GENERIC  OBJECT/BACKGROUND 
MODEL. 

The  general  modeling  strategy  described  in  Section  2  has  been  implemented  for  FLIR 
images  with  immediate  objectives  of  image  restoration  (i)  to  “smooth”  and  enhance 
homogeneous  subregions  corresponding,  for  example,  to  an  object  or  to  a  large  component 
of  an  object  of  interest,  and  (ii)  to  highlight  boundaries  between  separate  homogeneous 
subregions  as  a  precurser  to  object  detection  and  recognition.  We  have  designed  and 
implemented  a  two- level  hierarchical  MRF  model  combining  the  directly  observable  pixel 
process  and  a  linked  unobservable  binary  process  indicating  the  presence  or  absence  of 
edge  elements.  Models  like  the  one  described  here  were  suggested  and  illustrated  in  [2], 
3.1  Scene  Model. 

The  image  process  X  comprises  the  pixel  process  Xp  and  the  edge  process  Xs ,  X  = 
{.VP,X£}.  As  usual,  the  pixel  sites  form  an  R  x  C  array  of  points  (i?  rows  and  C 
columns)  in  a  square  lattice  arrangement.  We  denote  this  R  x  C  array  by  Sp .  The 
sites  for  the  edge  process,  collectively  denoted  SE ,  also  form  a  regular  graph  structure, 
envisioned  as  fitting  between  the  sites  in  Sp .  Let  u,v  denote  pixel  sites  in  the  square 
lattice  Sp,  For  each  pair  u,v  of  horizontally  or  vertically  adjacent  pixels,  there  exists 
an  “edge  site”  denoted  <  u,u  >  in  SE .  The  edge  site  s  =<  u,v  >  corresponds  to  the 
location  of  possible  edge  or  boundary  element  between  pixels  u  and  v.  The  edge  variables 
are  binary,  with  XEU  V>  equalling  1  or  0  to  indicate  the  presence  or  absence  of  an  edge 
element  at  <  u,v  >.  The  process  XE  consists  of  R[C  —  1)  +  C[R  —  1)  variables  XEUU>. 

The  totality  of  edge  and  pixel  sites  is  denoted  by  S.  (The  generic  point  s  may  refer  to 
a  pixel  or  to  an  edge  site  <  u,v  >.)  The  neighborhood  system  G  =  {G„s  6  5}  governs 
the  Markovian  dependence  structure  of  X  =  {Xp ,XE).  The  size  of  the  neighborhood 
determines  the  range  of  interactions.  We  restrict  our  design  of  the  process  to  “small” 
or  “local”  neighborhood  sets  G,,  to  keep  the  mathematical  models  as  simple  as  possible 
and  to  assure  feasibility  of  computational  procedures. 

We  adopt  the  following  neighborhood  system.  Each  pixel  site  has  eight  pixel  neigh¬ 
bors,  the  nearest  ones,  and  four  edge  neighbors.  Each  edge  site  <  u,v  >  has  six  edge 
neighbors — corresponding  to  possible  continuations  of  a  boundary  from  <  u,v  > — and 
the  two  pixel  neighbors  u  and  v.  Sites  near  the  boundary  of  the  lattice  have  fewer 
neighbors.  The  canonical  pixel  neighborhood  G,  and  edge  neighborhood  G<t  t>  are 
depicted  in  the  figure  below,  where  circles  represent  pixels  and  pluses  represent  edge 


sites.  (We  believe  this  edge  graph  originated  in  [11].) 
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To  illustrate  the  functional  form  of  the  models,  suppose  first  that  we  are  only  in¬ 
terested  in  modeling  “smoothness”  or  “regularity”  in  the  intensity  array  Xp ,  i.e.,  the 
tendency  of  nearby  pixels  to  have  similar  intensities.  Then  a  suitable  model  might  be 
X  =  Xp  with 

n(X  =  x)  =  Z~x  exp{0  J2  -  *«)}• 

(*.0 

where  the  sum  extends  over  all  neighboring  pairs  ( s,t )  of  pixels.  (Thus  each  interior 
pixel  is  included  in  eight  terms  in  the  summation.)  Here  4>  =  is  an  even  function, 
decreasing  for  6  >  0;  0  is  a  parameter  which  corresponds  to  “inverse  temperature”  and 
it  governs  the  degree  of  regularity.  It  is  distinct  from  the  “artificial  temperature”  T 
introduced  for  the  annealing  algorithm  (Section  2.5).  The  coefficient  is  introduced 
to  allow  different  weighting  of  pixel  pairs  oriented  in  different  directions.  We  commonly 
fix  C(i  {)  =  1  for  the  horizontal  and  vertical  pairs  and  =  l/\/2  for  diagonally 

adjacent  pairs.  A  renormalization  argument  shows  that  this  weighting  is  “asymptotically 
correct”  in  order  for  the  discrete  digital  images  Xp  to  approximate  rotationally  invariant 
(isotropic)  images  on  a  continuous  background  [8].  The  weights  also  permit  accurate 
modeling  of  anisotropic  FLIR  images. 

A  flexible  and  well-tested  choice  for  <f>  is  of  the  form 


(3.1) 


m  = 


where  B  is  a  parameter  depending  on  the  dynamic  range  of  the  image.  An  attractive 
feature  of  this  (^-function — in  contrast  to  one  that  decreases  without  bound — is  that  it 
does  not  attach  ever  increasing  penalties  to  larger  differences  6,  and  thus  it  will  allow 
sharp  gradients  in  intensity  across  region  boundaries.  A  choice  such  as  <£(<5)  =  —  62  would 
a  priori  inhibit,  indeed  prohibit,  adjacent,  internally  homogeneous  subregions  with  highly 
separated  intensities. 

With  the  inclusion  of  the  edge  process  XE  we  incorporate  our  expectations  about 
both  the  interactions  between  intensities  and  edges — i.e.,  where  the  edges  belong — and 
about  clusters  of  nearby  edges.  We  are  not  exactly  modeling  entire  boundaries  with  this 
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two-level  model,  but  rather  segments  of  boundaries;  except  in  the  simplest  imagery  and 
with  larger  neighborhoods,  it  is  essentially  impossible  to  distinguish  actual  boundary 
segments  from  intensity  gradients  due  to  lighting,  texture,  etc. 

For  the  pixel-edge  process,  the  complete  energy  function  U  =  U(XP,XE )  is  decom¬ 
posed  into  two  components: 

U{XP,XE)  =  U\XP,XE)  +  U\XE). 

We  construct  U 1  so  that  the  most  likely  configurations  will  have  XEt  t>  —  1  (respectively 
0)  when  the  intensity  difference  \xp  —  xp\  between  neighboring  pixels  is  large  (resp. 
small).  Put  differently,  we  want  to  “break”  the  bond  between  pixels  s  and  t  when  their 
values  are  “far”  apart.  Thus  we  choose 

(3.2)  U1(xp,xE)  =  -Y:0iC{,,t)(4>(xp  -xp)-02)  x  (1  -  IM(XE)) 

(».«) 

where  #i  >  &2  >  0.  The  value  of  6  for  which  C[a  t)<j>(6)  =  02  represents  an  intensity 
difference  for  which  we  have  “no  preference”  in  regard  to  the  on-off  state  of  an  edge;  such 
interpretations  of  the  model  parameters  are  helpful  when  one  is  setting  or  estimating 
values  of  the  parameters.  Finally,  in  equation  (3.2),  7(J  t)(X£)  =  1  when  the  XE  process 
“breaks”  the  bond  between  pixels  s  and  t,  and  /(4,t)(X£)  =  0  otherwise.  In  particular,  if 
s  and  t  are  horizontal  or  vertical  neighbors,  then  /(S)t)(X£)  =  XEa  t>  and  if  s  and  t  are 
diagonal  neighbors,  then  I^,  ^{XE)  is  a  Boolean  function  of  four  edge  elements  between 
s  and  t  requiring,  for  its  value  to  be  1,  that  at  least  two  of  the  edge  elements  are  “on  ' 
and  that  they  connect  to  form  a  segment  separating  s  from  t. 

The  remaining  component  U 2  of  the  total  energy  function  governs  the  organization 
of  nearby  edges.  We  define 

U2(xE)  =  -03J2Vd(xe) 

D 

where  63  >  0  and  where  the  sum  extends  over  all  subsets  D  of  four  neighboring  edge 
sites — the  maximal  “cliques”  in  the  edge  graph.  The  clique  function  Vp  assigns  weights 
in  accordance  with  our  expectations  about  edge  behavior.  More  specifically,  there  are 
six  possible  clique  states,  up  to  rotational  equivalence: 

0  |  °  00  o|o  o|o  00  00 

O  I  o  00  00  00  00  00 


Here  the  bars  indicate  that  the  edge  variable  at  the  indicated  site  is  “on”.  Let  VD  =  £,, 
for  t  =  1 , .  -  - ,  6,  denote  the  weights  assigned  to  the  six  configurations  above.  If  we 


assume  that  most  pixels  are  not  next  to  boundaries,  that  edges  should  continue,  and 
that  boundary  congestion  is  unlikely,  then  we  might  choose  £1  <  £2  <  $3  <  £<  <  $5  < 

A  specific  image-dependent  choice  is  made  in  the  experiment  described  in  Section  5. 

One  final  point  about  the  scene  model:  it  is  useful  to  write  the  total  energy,  up  to 
an  additive  constant,  as 

(3.3)  -U(x)  =  8^  C(,, )«*>(*?  -  xf)(l  -  /(,,«)(**))  +  d2  Y,  +  *3  £VD(xE) 

(*.0  (».«)  D 

For  inferential  purposes,  this  shows  that  our  model  is  an  exponential  family  in  8  = 
(01,02,03)-  In  addition,  the  form  in  (3.3)  is  useful  for  parameter  interpretation;  for 
instance,  it  becomes  clear  that  6-2  is  a  “reward”  for  edges. 

3.2  Degradations. 

The  Gibbs  distribution  determined  by  the  energy  function  U  in  equation  (3.3)  models 
the  ideal  scenes.  There  are  several  types  of  degradations  that  corrupt  an  ideal  scene 
before  it  is  observed.  Most  of  these  effects  are  well  understood  and  can  be  modeled 
accurately  in  terms  of  the  physical  processes  that  underlie  them.  In  the  end,  the  first 
approximation  of  the  degraded  observed  image  Y  will  reduce  to  the  pixel  process  Xp 
plus  additive  noise.  The  approximation  is  a  gross  simplification,  even  if  it  is  reasonably 
effective  as  a  basis  for  restoration  algorithms.  Ongoing  research  is  exploring  the  use  of 
more  accurate  degradation  models  which  incorporate  degradations  modeled  by  convo¬ 
lutions;  as  we  describe  below,  these  latter  degradations  include  atmospheric  absorption 
and  scattering,  diffraction  from  geometric  optics,  and  blurring  from  signal  averaging  and 
sampling  by  the  IR  sensor. 

Two  useful  references  for  understanding  degradations  of  IR  images  are  the  NV&EOL 
Technical  Report  [16]  by  J.A.  Ratches  et.  al.  and  the  NV&EOL  internal  working  paper 
[15]  prepared  for  our  project  by  V.  Mirelli.  Some  of  the  basic  physics  of  IR  radiation  and 
detection  is  described  in  [17]. 

The  primary  sources  of  IR  image  degradation  are: 

•  The  actual  thermal  radiation  from  the  ideal  scene  is  random  and  additive  to  Xp . 
The  random  component  has  mean  value  0  and  has  signal-dependent  variance  pro¬ 
portional  to  Xp .  The  exact  distribution  of  the  emitted  radiation  is  well-modeled 
by  a  Poisson  process  and  a  Gaussian  approximation  is  justified  by  convergence  of 
the  Poisson  Law  to  the  Normal. 


During  atmospheric  transmission  of  the  radiation,  there  is  absorption — dependent 
on  air  temperature  and  relative  humidity — and  scattering — dependent  on  visibility. 


The  scattering  is  normally  modeled  by  Beer’s  Law  [16].  The  effects  of  absorption 
and  scattering  enter  the  mathematical  model  in  the  form  of  a  convolution  of  the 
signal  with  a  kernel  depending  on  atmospheric  parameters  and  range-to-target. 

At  the  sensor,  the  first  degradation  stems  from  optical  diffraction.  The  geometrical 
optical  effect  is  modeled  by  a  convolution  of  the  signal  with  a  kernel  depending 
on  parameters  of  the  optical  system  (lens  diameter,  focal  lengths)  and  on  the 
wavelength  of  the  electromagnetic  radiation. 

Black-body  radiation  from  the  positive  temperature  of  the  detector  corrupts  the 
image  incident  at  the  detector.  This  effect  enters  the  model  as  additive  “noise”  on 
top  of  the  signal. 

The  electromagnetic  energy  in  the  IR  radiation  is  converted  to  an  electrical  response 
by  the  sensor.  The  response  is  a  random  process  subordinated  on  the  input.  This 
can  be  represented  mathematically  as  signal-dependent  additive  noise,  again  with  a 
Poisson  distribution,  where  both  the  conditional  mean  and  variance  of  the  response 
(given  the  input)  equals  the  input. 

The  electrical  response  is  digitized  through  a  combination  of  averaging  and  sam¬ 
pling.  Conceptually,  a  scanning  detector  returns  a  continuous  response  which  is 
averaged  in  the  direction  of  the  scan  and  which  is  discretely  sampled  in  the  direction 
orthogonal  to  the  scan.  The  combination  of  averaging  and  sampling  implies  that 
the  observed  process  will  not  be  isotropic.  Digitization  is  described  mathematically 
through  a  convolution  of  the  continuous  signal  with  a  singular  kernel. 

Finally,  electronic  noise  may  enter  at  the  last  stage  of  actually  observing  the  digi¬ 
tized  signal.  The  noise  enters  as  an  additive  effect,  independent  of  the  signal. 


4  IMPLEMENTATION  OF  THE  RESTORATION 
ALGORITHM. 


The  following  subsections  give  a  complete  listing  of  a  standard  FORTRAN77  program 
that  implements  stochastic  relaxation,  with  optional  annealing,  for  the  model  described 
in  Section  3.  The  subroutine  that  computes  the  dependence  of  the  total  energy  on  the 
edge  process  (SUBROUTINE  DEE)  actually  implements  a  model  that  is  slightly  more 
general  than  equation  (3.3).  It  incorporates  a  parameter  (CE2C)  which  inhibits  the 
occurrence  of  nearby  parallel  edges.  The  model  of  Section  3  is  implemented  by  this 
program  when  CE2C=0. 

This  program  has  been  delivered  to  the  Advanced  Modeling  Team  at  NV&EOL  and 
has  been  used  there  for  experiments  with  restoration  of  FLIR  images.  The  presumptions 
about  formats  of  input  and  output  files  are  best  documented  by  the  input  and  output 
subroutines  READIN  and  WRITEO,  which  are  listed  below.  Experiments  with  the  use 
of  this  program  are  described  in  Section  5. 

4.1  Main  Program  RESTOR. 

The  main  program  guides  input,  output  and  stochastic  relaxation  of  the  pixel  and  edge 
processes. 


PROGRAM  RESTOR 
C  SET  UP  DATA  STRUCTURES 
INCLUDE  (COMMON) 

C  TYPES 

INTEGER  NIT 
C  GET  INPUT 

CALL  READIN 
C  ITERATE 

DO  10  NIT-NSTART.NSTOP 

PRINT  *.  ’ITERATION  NIT 
IF  (NIT.LE.NO)  THEN 
T-TO 
ELSE 

T*T0/(1 . 0+L0G(FL0AT(NIT-N0))) 
END  IF 

PRINT  *.  ’TEMPERATURE  ’.  T 
IF  (IXP.Eq.l)  CALL  ITXP 
IF  (IXE.EQ.l)  CALL  ITXE 
10  CONTINUE 
C  OUTPUT  RESULTS 


RES00010 

RES00020 

RES00030 

RES00040 

RES00050 

RES00060 

RES00070 

RES00080 

RES00090 

RES00100 

RES00110 

RES00120 

RES00130 

RES00140 

RES00150 

RES00160 

RES00170 

RES00180 

RES00190 

RES00200 


4.2  Include  File  COMMON. 


The  “include”  file  declares  global  variables,  sets  parameter  values,  and  defines  COMMON 
blocks. 

CE1A  is  the  model  parameter  0X  (equation  3.3). 

CElB  is  the  model  parameter  B  (equation  3.1). 

CE2A  is  the  model  parameter  03  (equation  3.3). 

CE2B  is  the  model  parameter  (equation  3.3). 

CE2C  is  not  used  in  model  (3.3)  and  is  set  to  0. 

PMIN  is  the  minimum  value  of  the  range  of  the  pixel  process  if. 

PMAX  is  the  maximum  value  of  the  range  of  the  pixel  process  if. 

SIGMA  is  the  standard  deviation  of  the  additive  noise  corrupting  the  observed 
process  F. 

MAXDA  is  the  maximum  number  of  equally  spaced  discrete  levels  used  to  quantize 
the  range  [PMIN, PM  AX]  of  if . 

NDA  is  the  actual  number  of  equally  spaced  discrete  levels  used  to  quantize  the  range 
[PMIN, PMAX]  of  if. 


C  CONSTANTS 

INTEGER  NX,  NY.  MAXDA 
REAL  DIAG 

PARAMETER  (NX-64 .NY-84 .MAXDA-100, DIAG- . 707) 

C  DECLARE  PARAMETERS.  WHICH  WILL  BE  READ  FROM  UNIT  7 
REAL  CE1A . CElB . CE2A .CE2B . CE2C .PMIN , PMAX . SIGMA 
C  VARIABLES  AND  ARRAYS 

INTEGER  IS.  ID,  IP,  IXP,  IXE.  NO.  NSTART,  NSTOP. 

M  NDA 

REAL  T0,T,XP(0:NX+1 ,0:NY+1) . XE(-1 : NX+2 , -1 : NY+2 .2) . XP0(NX .NY) , 
M  ADSIG , SIGSQD 

DOUBLEPRECISION  SEED 
C  COMMON  GLOBAL  DATA  STRUCTURES 

COMMON  SEED . CE1A . CElB , CE2A , CE2B , CE2C . PMIN , PMAX .SIGMA . 

M  TO. T.XP.XE.XPO, ADSIG. SIGSqD, 

M  IS.ID.IP.IXP.IXE.NO. NSTART . NSTOP . NDA 


COMOOOIO 
C0MOOO2O 
C0MOOO3O 
CDM00040 
COMOOOBO 
COM00060 
CQM00070 
COM00080 
C0M00090 
COMOOIOO 
C0M00110 
COMOO 120 
COMOO 130 
C0M00140 
COMOO 150 
COMOO 160 


4.3  Subroutine  READIN 


The  input  routine  READIN  prompts  the  user  for  interactive  input  of  program  and  model 
parameters  and  reads  in  files  containing  images,  including  the  observed  image  and  any 
results  that  may  be  available  from  previous  processing  by  the  relaxation  algorithm. 


SUBROUTINE  READIN 
C  SET  UP  DATA  STRUCTURES 
INCLUDE  (COMMON) 

C  TYPES 

INTEGER  I,  J.  K 
C  EXTERNAL  FUNCTIONS  CALLED 
REAL  GGUBFS 

C  READ  PARAMETER  VALUES  FROM  UNIT  7 
READ (7,*).  CE1A 
READ (7,*).  CE1B 
READ (7 . *)  .  CE2A 
READ (7,*),  CE2B 
READ(7 , *) ,  CE2C 
READ (7.*),  PMIN 
READ (7.*),  PMAX 
READ (7.*),  SIGMA 
CLOSE (UN IT-7) 

C  DETERMINE  IF  GOAL  IS  IMAGE  SAMPLING 
IS-0 

WRITE(6 . *) ,  ’ENTER  1  IF  A  SAMPLE  IMAGE  IS  DESIRED,  0  IF  PURPOSE’ 
WRITE(0 , *) .  ’IS  RESTORATION’ 

READ (6,*),  IS 

C  IF  GOAL  IS  RESTORATION.  DETERMINE  IF  ORIGINAL  IMAGE  IS  RESULT  OF 
C  A  DEGRADATION 
ID-0 

IF  (IS.Eq.O)  THEN 

WRITE(0 .  *)  ,  ’ENTER  1  IF  THERE  IS  A  DEGREDATION,  0  OTHERWISE’ 
READ (5 . *) ,  ID 
END  IF 

C  DETERMINE  IF  IMAGE  HAS  ALREADY  BEEN  PARTIALLY  PROCESSED 
IP-0 

WRITE(@, *) ,  ’ENTER  1  IF  PROCESSING  BEGAN  WITH  A  PREVIOUS  RUN,’ 
WRITE(6 , *) ,  ’0  OTHERWISE’ 

READ(6 , *) ,  IP 

C  DETERMINE  WHICH  LEVELS  OF  HIERARCHY  ARE  TO  BE  ACTIVE 
IXP-0 


REA00010 

REA00020 

REA00030 

REA00040 

REA00050 

REA00060 

REA00070 

REA00080 

REA00090 

REA00100 

REA00110 

REA00120 

REA00130 

REA00140 

REA00150 

REA00160 

REA00170 

REA00180 

REA00190 

REA00200 

REA00210 

REA00220 

REA00230 

REA00240 

REA00250 

REA00260 

REA00270 

REA00280 

REA00290 

REA00300 

REA00310 

REA00320 

REA00330 

REA00340 

REA00350 

REA00380 
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IXE=0  REA00370 

WRITE(8 , *) ,  'ENTER  1  IF  PIXEL  PROCESS  WILL  BE  ACTIVE,  0  OTHERWISE 'REA00380 
READ (5,*),  IXP  REA00390 

WRITE(6 ,  *)  .  ’ENTER  1  IF  EDGE  PROCESS  WILL  BE  ACTIVE,  0  OTHERWISE'  REA00400 
READ (6,*)  IXE  REA00410 

C  DETERMINE  NUMBER  OF  DISCRETE  VALUES  REA00420 

WRITE(6 , *) ,  'ENTER  NUMBER  OF  GREY  LEVELS’  REA00430 

WRITE(6 , *) .  '(NO  MORE  THAN '.MAXDA,')'  REA00440 

READ(5 , *)  .  NDA  REA00450 

C  DETERMINE  TEMPERATURE  CONTROL  PARAMETERS  REA00460 

WRITE(6 , *) ,  'ENTER  STARTING  TEMPERATURE,  EVEN  IF'  REA00470 

WRITE(6 , *)  .  'THIS  IS  FROM  A  PREVIOUS  RUN’  REA00480 

RE AD (5,*),  TO  REA00490 

WRITE(6 , *)  ,  'ENTER  NUMBER  OF  ITERATIONS  BEFORE  INITIATION'  REA00500 

WRITE(6 ,  *)  ,  ’OF  ANNEALING'  REA00510 

READ (5 , *) ,  NO  REA00520 

C  DETERMINE  STARTING  AND  STOPPING  ITERATIONS  REA00630 

WRITE (6 , *) ,  'ENTER  NUMBER  OF  FIRST  ITERATION  FOR  THIS  RUN’  REA00B40 

READ ( 5 . * )  ,  NSTART  REA00660 

WRITEC6,*),  'ENTER  NUMBER  OF  LAST  ITERATION  FOR  THIS  RUN’  REA00660 

READ(5 . *) .  NSTOP  REA00570 

C  GET  SEED  FOR  RANDOM  NUMBER  GENERATOR  REA00B80 

WRITE(6 , *) ,  'ENTER  SEED  FOR  RANDOM  NUMBER  GENERATOR’  REA00B90 

RE AD ( B . * )  .  SEED  REA00600 

C  IF  GOAL  IS  RESTORATION.  AND  THERE  IS  A  DEGRADATION,  THEN  REA00610 

C  DETERMINE  THE  STANDARD  ERROR  OF  ANY  NOISE  THAT  HAS  BEEN  ADDED  TO  REA00620 

C  THE  IMAGE  AND  COMPUTE  THE  TOTAL  SIGMA  SQUARED  ("SIGSQD")  REA00630 

IF  (IS.Eq.O.AND.ID.Eq.l)  THEN  REA00640 

WRITE(6 , *) ,  'ENTER  STANDARD  ERROR  OF  ADDED  NOISE  (0  IF  NO’  REA006B0 

WRITE(6 , *) ,  'NOISE  HAS  BEEN  INTRODUCED)'  REA00660 

READ(6 , *) ,  ADSIG  REA00670 

SIGSQD«ADSIG**2+SIGMA**2  REA00680 

END IF  REA00690 

C  READ  IN  DATA  REA00700 

IF  (IP.Eq.l)  THEN  REA00710 

DO  1  J-l.NY  REA00720 

READ(l.e)  (XP(I.J).I-l.NX)  REA00730 

1  CONTINUE  REA00740 

DO  3  K*1 ,2  REA00760 

DO  4  J-l.NY  REA00760 

READU.6)  (XE(I.J.K)  .1-1, NX)  REA00770 


-*V  *- 
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4  CONTINUE 

3  CONTINUE 

6  FORMAT (10F7. 2) 

CLOSE(UNIT-l) 

END  IF 

IF  (ID.Eq.l)  THEN 
DO  7  J-l.NY 

READ(2 ,6)  (XPO(I.J).I-l.NX) 

7  CONTINUE 
CL0SE(UNIT*2) 

ENDIF 

IF  (IS.EQ.O.AND. ID.EQ.O. AND. IP.EQ.O)  THEN 
DO  0  J-l.NY 

READ (3 ,6)  (XP(I.J) .1-1. NX) 

9  CONTINUE 

CL0SE(UNIT-3) 

ENDIF 

C  INITIALIZE  DATA  ARRAYS.  ALL  NONPIXEL  STRUCTURES  ARE 
C  INITIALIZED  TO  "NOT  PRESENT".  UNLESS  THERE  WAS 
C  PREVIOUS  PROCESSING. 

IF  ( ID . EQ . 1 . AND .IP.EQ.O)  THEN 
DO  16  1*1 .NX 
DO  20  J-l.NY 

XP(I.J)-XPO(I.J) 

20  CONTINUE 

15  CONTINUE 

ENDIF 

IF  ( IS . EQ . 1 . AND  IP.EQ.O)  THEN 
DO  60  1*1. NX 
DO  70  J*1 , NY 

XP ( I . J ) -PMIN+ (PMAX - PMIN ) *GGUBFS ( SEED ) 

70  CONTINUE 

60  CONTINUE 

ENDIF 

IF  (IP.EQ.O)  THEN 
DO  76  K-1.2 
DO  80  J-l.NY 
DO  90  I-l.NX 
XE(I . J , K) *0 . 0 
90  CONTINUE 

80  CONTINUE 


I'VvvvV 


REA00780 
REA00790 
REA00800 
REA00810 
REA00820 
REA00830 
REA00840 
REA00850 
REA00860 
REA00870 
REA00880 
REA00890 
REA00900 
REA00910 
REA00920 
REA00930 
REA00940 
REA 00960 
REA00960 
REA00970 
REA00980 
REA00990 
REA01000 
REA01010 
REA01020 
REA01030 
REA01040 
REA01050 
REA01060 
REA01070 
REA01080 
REA01090 
REA01100 
REA01110 
REA01120 
REA01130 
REA01140 
REA01150 
REA01160 
REA01170 
REA01180 
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75  CONTINUE 

ENDIF 

C  INITIALIZE  DUMMY  BOUNDARIES 
DO  100  J=0 , NY+1 
XP(0,J) =1000.0 
XP (NX+1 , J)=1000 . 0 
100  CONTINUE 

DO  110  1*1 , NX 
XP(I.0)=1000.0 
XP(I , NY+1 )=1000 . 0 
110  CONTINUE 

DO  120  I=-l ,NX+2 
XE(I.-l,l)-0.0 
XE(I,-1 ,2) =0.0 
XE  ( 1 , 0 , 1 )  =0 . 0 
XE(I ,0 , 2) =0 . 0 
XE(I ,NY , 2)=0 . 0 
XE(I , NY+1 , 1) =0 . 0 
XE(I , NY+1 ,2) =0 . 0 
XE( I ,NY+2 , 1) =0 . 0 
XE(I ,NY+2 .2) =0 . 0 
120  CONTINUE 

DO  130  J--1.NY+2 
XE(-1,J,1)=0.0 
XE(- 1 , J , 2)=0 . 0 
XE(0 , J , 1) =0 . 0 
XE(0 , J , 2) =0 . 0 
XE(NX , J , 2)=0 . 0 
XE(NX+1 , J  ,  1 ) =0 . 0 
XE(NX+1 , J ,2) =0 . 0 
XECNX+2, J,1)=0.0 
XE(NX+2 . J ,2) =0 . 0 
130  CONTINUE 
END 


REA01190 

REA01200 

REA01210 

REA01220 

REA01230 

REA01240 

REA01250 

REA01260 

REA01270 

REA01280 

REA01290 

REA01300 

REA01310 

REA01320 

REA01330 

REA01340 

REA01350 

REA01360 

REA01370 

REA01380 

REA01390 

REA01400 

REA01410 

REA01420 

REA01430 

REA01440 

REA01450 

REA01460 

REA01470 

REA01480 

REA01490 

REA01500 

REA01510 

REA01520 


4.4  Subroutine  WRITEO. 

The  output  routine  WRITEO  writes  the  output  image  file  to  the  disk. 


SUBROUTINE  WRITEO 

WRIOOOIO 

c 

SET  UP  DATA  STRUCTURES 

WRI00020 

INCLUDE  (COMMON) 

WRI00030 

c 

TYPES 

WRI00040 

INTEGER  I,  J,  K 

WRIOOOBO 

c 

WRITE  OUTPUT  TO  UNIT  4 

WR I 00060 

DO  1  J-l.NY 

WRI00070 

WRITE(4.6)  (XP(I.J) .1-1. NX) 

WR I 00080 

1 

CONTINUE 

WR I 00090 

DO  3  K-1,2 

WRI00100 

DO  4  J-l.NY 

WRI001 10 

WRITE(4 ,6)  (XE(I.J.K) ,I«1.NX) 

WRI00120 

4 

CONTINUE 

WRI00130 

3 

CONTINUE 

WRI00140 

6 

FORMAT (10F7. 2) 

WRI00160 

CL0SE(UNIT«4) 

WRI00160 

END 

WRI00170 

4.5  Subroutine  ITXP 


The  subroutine  ITXP  guides  the  execution  of  the  relaxation  algorithm  for  the  pixel 
process  Xp . 


SUBROUTINE  ITXP 
C  SET  UP  DATA  STRUCTURES 
INCLUDE  (COMMON) 

C  TYPES 

INTEGER  I.J.K 

REAL  EP(MAXDA) ,  SUM(MAXDA) .  TOT,  EMIN.  EMAX .  NRAND 
C  EXTERNAL  FUNCTIONS  CALLED 
REAL  GGUBFS 
C  ITERATE  PIXEL  VALUES 
DO  10  J-l.NY 
DO  20  1*1 . NX 

C  COMPUTE  ENERGY  VECTOR  FOR  PIXEL  (I.J)  AND  STORE  IN  EP  EP(K) 

C  IS  THE  RELATIVE  ENERGY  FOR  XP(I.J)  AT  THE  K'TH  DISCRETE  VALUE 
CALL  PIXEN(I.J.EP) 

C  PREVENT  OVERFLOWS  AND  UNDERFLOWS  BY  RESCALING  AND  TRUNCATING  EP 
EMIN  *EP(1) 

DO  5  K-2.NDA 

IF  (EP(K)  LT.EMIN)  THEN 
EMIN-EP(K) 

END  IF 

5  CONTINUE 
EMAX*T»20 . 0 
DO  6  K-l.NDA 

EP (K) *MIN(EMAX ,EP(K) -EMIN) 

6  CONTINUE 

C  UPDATE  VALUE  OF  XP(I.J) 

SUM(1)*EXP(-EP(1)/T) 

DO  30  K-2.NDA 

SUM(K)*SUM(K-l)fEXP(-EP(K)/T) 

30  CONTINUE 

NRAND*GGUBFS (SEED ) 

TOT-SUM (ND A) 

DO  40  K-l.NDA 

IF  (NRAND. LE  (SUM(K) /TOT))  THEN 

XP( I , J)*PMIN* ( (PMAX-PMIN) * (K-1))/(NDA-1) 

GO  TO  20 

ENDIF 


ITX00010 
ITX00020 
ITX00030 
ITX00040 
ITX00060 
ITX00060 
ITX00070 
ITX00080 
ITXOOOOO 
ITX00100 
ITX001 10 
ITX001 20 
ITX00130 
ITX00140 
ITXC0160 
ITX00160 
ITX00170 
ITX00180 
ITX00100 
ITX00200 
ITX00210 
ITX00220 
ITX00230 
ITX00240 
ITX00250 
IT 100260 
ITI00270 
ITX00280 
ITX00200 
ITX00300 
ITX00310 
ITI00320 
ITX00330 
ITX00340 
ITX00350 
I T 1 00 360 
ITX00370 
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40 

20 

10 


CONTINUE 

CONTINUE 

CONTINUE 

END 


ITX00380 

ITX00390 

ITX00400 

ITX00410 


4.6  Subroutine  PIXEN 


The  subroutine  PIXEN  is  called  by  ITXP  and  returns  the  vector  of  (relative)  energies 
that  determine  the  local  conditional  distribution  of  the  possible  values  for  the  pixel 
process  at  an  arbitrary  site. 


C  PIXEN(I.J.EP)  COMPUTES  THE  RELATIVE  ENERGY  FOR  THE  NDA  DIFFERENT 
C  POSSIBLE  LEVELS  OF  PIXEL  (I.J).  THESE  ARE  RETURNED  VIA  EP(MAXDA) . 

SUBROUTINE  PIXENCI.J.EP) 

C  SET  UP  DATA  STRUCTURES 
INCLUDE  (COMMON) 

C  TYPES 

INTEGER  I,  J,  K 

REAL  EP(MAXDA) ,  ADIFF,  XPTEMP ,  INC 
C  INITIALIZE  EP 

DO  10  K-l.NDA 
EP(K)=0.0 
10  CONTINUE 

C  COMPUTE  DEGRADATION  CONTRIBUTION  TO  ENERGY  (IF  ANY) 

IF  (ID.Eq.l)  THEN 

CALL  PIXENO(I.J.EP) 

END  IF 

C  COMPUTE  PURE  PIXEL  CONTRIBUTION  TO  ENERGY 
INC- (PMAX-PMIN) / (NDA- 1 ) 

DO  20  K-l.NDA 

XPTEMP-PMIN+INC* (K- 1 ) 

C  PIXEL  TO  UPPER  LEFT: 

IF  ((XECI-1.J, 1) +XE(I-1 , J-l , 2) ) * 

M  (XE(I-1 , J-l , 1)+XE(I . J-1,2) ) .LT. . 5)  THEN 
ADIFF-ABS((XPTEMP-XP(I-1, J-1))/CE1B) 

EP(K)-EP(K) -CE1A*DIAG/ (1 . 0+ ADIFF* AD IFF) 

ENDIF 

C  PIXEL  ABOVE: 

IF  (XE(I , J-l .2) .LE . .5)  THEN 

ADIFF»ABS((XPTEMP-XP(I , J- 1) )/CElB) 
EP(K)-EP(K)-CE1A/(1.0+ADIFF*ADIFF) 

ENDIF 

C  PIXEL  TO  UPPER  RIGHT: 

IF  ( (XE(I , J-l ,2)+XE(I . J-l , 1) ) * 

M  (XE(I , J . 1)+XE(I+1 , J- 1 ,2) ) .LT . .5)  THEN 
ADIFF-ABS( (XPTEMP -XP(I+1 , J- 1) ) /CE1B) 

EP(K) -EP(K) -CE1A*DIAG/ ( 1 . 0+ AD IFF* ADIFF) 


PIX00010 
PIX00020 
PIX00030 
P 1X00040 
P 1X00050 
P 1X00060 
P 1X00070 
P 1X00080 
P 1X00090 
PIX00100 
PIX00110 
PIX00120 
PIX00130 
P 1X00140 
P 1X00150 
PIX00160 
P 1X00170 
PIX00180 
PIX00190 
PIX00200 
PIX00210 
PIX00220 
PIX00230 
PIX00240 
P 1X00250 
P 1X00260 
PIX00270 
P 1X00280 
P 1X00290 
P 1X00300 
PIX00310 
P 1X00320 
P 1X00330 
P 1X00340 
P 1X00350 
P 1X00360 
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END  IF 

C  PIXEL  TO  LEFT: 

IF  (XE(I-l.J.l) .LE. .5)  THEN 

ADIFF-ABS((XPTEMP-XP(I-1.J))/CE1B) 

EP (K) -EP (K) -CE1A/ ( 1 . 0+ADIFF*ADIFF) 

END  IF 

C  PIXEL  TO  RIGHT: 

IF  (XE(I.J.l).LE. .5)  THEN 

ADIFF-ABS((XPTEMP-XP(I+1,J))/CE1B) 
EP(K)-EP(K) -CE1A/(1 . 0+ADIFF*ADIFF) 

END  IF 

C  PIXEL  TO  LOWER  LEFT: 

IF  ( (XE(I-1 , J,2)+XE(I-1 , J , 1) ) * 

M  (XE(I-1 , J+l , 1)+XE(I . J ,2) ) .LT . .5)  THEN 
ADIFF-ABS((XPTEMP-XP(I-1,J+1))/CE1B) 
EP(K)-EP(K) -CE1A*DIAG/(1 . 0+ADIFF*ADIFF) 
END  IF 

C  PIXEL  BELOW: 

IF  (XE(I.J,2).LE. .5)  THEN 

ADIFF«ABS((XPTEMP-XP(I,J+1))/CE1B) 

EP (K) -EP (K) -CE1 A/ ( 1 . 0+ADIFF*ADIFF) 

END  IF 

C  PIXEL  TO  LOWER  RIGHT: 

IF  ( (XE(I , J , 2)+XE(I , J+l , 1) )* 

M  (XE(I,J.1)+XE(I+1.J.2)).LT. .5)  THEN 
ADIFF-ABS((XPTEMP-XP(I+1 , J+l) )/CElB) 
EP(K)-EP(K)-CE1A*DIAG/(1.0+ADIFF*ADIFF) 
ENDIF 

20  CONTINUE 
END 


PIX00370 
PIX00380 
P 1X00390 
PIX00400 
PIX00410 
P 1X004 20 
P 1X004 30 
P 1X00440 
PIX00460 
PIX00460 
P 1X00470 
PIX00480 
P 1X00490 
PIX00500 
PIX00510 
PIX00620 
PIX00530 
P 1X00540 
P 1X00550 
P 1X00560 
P 1X00570 
P 1X00580 
P 1X00590 
P 1X00600 
PIX00610 
P 1X00620 
P 1X00630 
PIX00640 
PIX00650 
P 1X00660 


4.7  Subroutine  PIXENO 


The  subroutine  PIXENO  is  called  by  PIXEN  and  returns  that  part  of  the  local  energy 
vector  attributable  to  the  difference  between  the  observed  image  and  the  current  state 
of  the  restoration. 


C  PIXENOU .J.EP)  COMPUTES  THE  DEGRADATION  CONTRIBUTION  TO 
C  THE  RELATIVE  ENERGY  FOR  THE  NDA  DIFFERENT  POSSIBLE  LEVELS 
C  OF  PIXEL  (I.J).  THESE  ARE  RETURNED  VIA  EP(MAXDA) 
SUBROUTINE  PIXENO( I , J . EP) 

C  SET  UP  DATA  STRUCTURES 
INCLUDE  (COMMON) 

C  TYPES 

INTEGER  I.  J.  K 

REAL  EP(MAXDA) .  XPTEMP .  INC.  TSIGSQ 
C  COMPUTE  DEGREDATION  CONTRIBUTION  TO  ENERGY 
INC*(PMAX-PMIN) / (NDA- 1 ) 

TSIGSQ«2*SIGSqD 
DO  10  K-l.NDA 

XPTEMP*PMIN-»INC*  (K- 1 ) 

EP(K) «EP(K) ♦ (XPTEMP -XPO( I . J) ) **2/TSIGSq 
10  CONTINUE 
END 


PIX00010 
P 1X00020 
P 1X00030 
P 1X00040 
PIX00050 
P 1X00060 
P 1X00070 
P 1X00080 
P 1X00090 
PIX00100 
PIXOOUO 
PIX00120 
PIX00130 
PIX00140 
PIX00150 
PIX00160 
PIX00170 


4.8  Subroutine  ITXE 


The  subroutine  ITXE  guides  the  execution  of  the  relaxation  algorithm  for  the  edge 
process  A'£. 


SUBROUTINE  ITXE 

ITX00010 

C  SET  UP  DATA  STRUCTURES 

ITX00020 

INCLUDE  (COMMON) 

ITX00030 

C  TYPES 

ITX00040 

INTEGER  I.  J.  K 

ITX00050 

REAL  PON. EXPO 

ITX00060 

C  EXTERNAL  FUNCTIONS  CALLED 

ITX00070 

REAL  DEE 

ITX00080 

C  ITERATE  EDGE  PROCESS 

ITX00090 

DO  10  K-l  .2 

ITX00100 

DO  20  J-l.NY+l-K 

ITX001 10 

DO  30  1*1 , NX-2+K 

ITX00120 

EXPO-MIN(10 . 0 .MAX ( - 10.0 ,DEE( I . J , K) /T) ) 

ITX00130 

PON-1/ (1+EXP (EXPO)) 

ITX00140 

IF  (GGUBFS(SEED)  LE.PON)  THEN 

ITX00150 

XE( I , J , K)-l  0 

ITX00160 

ELSE 

ITX00170 

XE(I , J ,K)-0 .0 

ITX00180 

END  IF 

ITX00190 

30  CONTINUE 

ITX00200 

20  CONTINUE 

ITX00210 

10  CONTINUE 

ITX00220 

END 

ITX00230 

4.9  Subroutine  DEE 


The  subroutine  DEE  is  called  by  ITXE  and  computes  the  energy  difference  between  the 
states  “on”  and  “off”  for  the  edge  element  at  an  arbitrary  edge  site. 


C  DEE  CALCULATES  THE  ENERGY  DIFFERENCE  (DELTA  ENERGY)  BETWEEN 
C  EDGE  ELEMENT  (I.J.K)  IN  STATE  1  (ON)  AND  EDGE  ELEMENT  (I.J.K) 
C  IN  STATE  0  (OFF)  . 

REAL  FUNCTION  DEE( I.J.K) 

C  SET  UP  DATA  STRUCTURES 
INCLUDE  (COMMON) 

C  TYPES 

INTEGER  I.  J.  K.  NON 
REAL  HOLD.  RON.  AD IFF 
C  INITIALIZE  DEE 
DEE-0.0 

C  COMPUTE  NONDIAGONAL  PIXEL/EDGE  CONTRIBUTION 
ADIFF*ABS( (XP(I.J)-XP(I*2-K,J+K-1))/CE1B) 

DEE-DEE+CE1A/(1  0+ADIFF*ADIFF) 

C  COMPUTE  NONDIAGONAL  "BOND-BREAKING"  PENALTY 
DEE-DEE-CE2B 

C  COMPUTE  4-CLIQUE  TERMS.  INCLUDING  DIAGONAL  PIXEL/EDGE 
C  TERMS  AND  DIAGONAL  BOND-BREAKING  TERMS 
H0LD-XE( I.J.K) 

XE ( I . J , K) -1 . 0 

IF  (K.Eq.l.AND. J.GT.l)  THEN 

R0N-XE(I , J-l , 1)+XE(I+1 , J-l ,2)+XE(I , J , 1)+XE(I , J-l ,2) 
NON-NINT(RON) 

IF  (NON.Eq.l)  THEN 
DEE-DEE+3*CE2A 
ELSEIF  (NON . EQ . 2)  THEN 
DEE-DEE-2*CE2A 
IF  (XE(I, J-l .2)  GT. .5)  THEN 
DEE-DEE-CE2B 

ADIFF-ABS( (XP( I ,J)-XP(I+1,J-1))/CE1B) 
DEE-DEE+CE1A*DIAG/ ( 1 . 0+ADIFF*ADIFF) 

ELSEIF  (XE(I+1 . J-l .2) .GT. .5)  THEN 
DEE-DEE-CE2B 

ADIFF-ABS( (XP(I ,J-1)-XP(I+1,J))/CE1B) 
DEE-DEE+CE1A*DIAG/ ( 1 . O+ADIFF* ADIFF) 

ELSE 

DEE«DEE-2*CE2B 


DEEOOOIO 

DEE00020 

DEE00030 

DEE00040 

DEE00050 

DEE00060 

DEE00070 

DEE00080 

DEE00090 

DEEOOIOO 

DEEOOllO 

DEE00120 

DEE00130 

DEE00140 

DEE00150 

DEE00160 

DEE00170 

DEE00180 

DEE00190 

DEE00200 

DEE00210 

DEE00220 

DEE00230 

DEE00240 

DEE00250 

DEE00260 

DEE00270 

DEE00280 

DEE00290 

DEE00300 

DEE00310 

DEE00320 

DEE00330 

DEE00340 

DEE00360 

DEE00360 

DEE00370 
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ADIFF-ABS( (XP(I ,J)-XP(I+ltJ-l))/CElB) 
DEE-DEE+CE1A*DIAG/(1.0+ADIFF*ADIFF) 
ADIFF-ABS( (XP(I ,J-1)-XP(I+1.J))/CE1B) 
DEE-DEE+CE1A*DIAG/ ( 1 . 0+ADIFF*ADIFF) 

END  IF 

ELSEIF  (N0N.EQ.3)  THEN 
DEE-DEE+CE2A 

IF  (XEd+l.J-1.2)  .LT.  .5)  THEN 
DEE-DEE-CE2B 

ADIFF-ABS((XP(I,J)-XP(I+1,J-1))/CE1B) 
DEE-DEE+CE1A*DIAG/ ( 1 . 0+ADIFF*ADIFF) 

ELSEIF  (XECI.J-1.2)  LT. .6)  THEN 
DEE-DEE-CE2B 

ADIFF-ABS((XP(I.J-1)-XP(I+1,J))/CE1B) 
DEE-DEE+CE1A*DI AG/ ( 1 . O+ADIFF* ADIFF) 

ENDIF 

ELSEIF  (NON.Eq.4)  THEN 
DEE-DEE+CE2A 
ENDIF 
ENDIF 

IF  (K . Eq . 1 . AND . J . LT . NY)  THEN 

RON*XE(I . J , 1) +XE(I+1 , J . 2) +XE(I ,  J*1 , 1)+XE(I , J ,2) 
NON-NINT(RON) 

IF  (NON.Eq.l)  THEN 
DEE-DEE+3*CE2A 
ELSEIF  (NON.Eq.2)  THEN 
DEE-DEE- 2*CE2A 
IF  (XE(I , J ,2) .GT. .5)  THEN 
DEE-DEE -CE2B 

ADIFF-ABS((XP(I. J)-XP(I+1.J*1))/CE1B) 
DEE-DEE+CE1A*DIAG/(1 . 0+ADIFF*ADIFF) 
ELSEIF  (XE(I+1.J,2) .GT. .6)  THEN 
DEE-DEE -CE2B 

ADIFF-ABS( (XP(I , J+l) -XP(I*1 , J) )/CElB) 
DEE-DEE+CE1 A*D I AG/ ( 1 . O+AD IFF  *ADIFF ) 

ELSE 

DEE-DEE-2*CE2B 

ADIFF-ABS((XP(I , J)-XP(I+1 ,J+1))/CE1B) 
DEE-DEE+CE1A*DIAG/(1 . 0+ADIFF*ADIFF) 
ADIFF-ABS( (XP(I , J+l) -XP(I+1 , J) )/CElB) 
DEE-DEE+CE1A*DIAG/ ( 1 . 0+ADIFF*ADIFF) 


DEE00380 

DEE00390 

DEE00400 

DEE00410 

DEE00420 

DEE00430 

DEE00440 

DEE00450 

DEE00460 

DEE00470 

DEE00480 

DEE00490 

DEE00500 

DEE00510 

DEE00520 

DEE00530 

DEE00540 

DEE00550 

DEE00580 

DEE00B70 

DEE00S80 

DEE00690 

DEE00600 

DEE00610 

DEE00620 

DEE00830 

DEE00640 

DEE00650 

DEE00660 

DEE00670 

DEE00680 

DEE00690 

DEE00700 

DEE00710 

DEE00720 

DEE00730 

DEE00740 

DEE00750 

DEE00760 

DEE00770 

DEE00780 


■■ 


END  IF 

ELSEIF  (N0N.EQ.3)  THEN 
DEE-DEE+CE2A 

IF  (XE(I+1 , J , 2) .LT. .5)  THEN 
DEE*DEE-CE2B 

ADIFF*ABS( (XP(I ,J)-XP(I+1,J+1))/CE1B) 
DEE»DEE+CE1A*DIAG/(1 . 0+ADIFF*ADIFF) 

ELSEIF  (XEd.J.2)  .LT.  .5)  THEN 
DEE-DEE-CE2B 

ADIFF-ABS((XP(I,J+l)-XPa+l.J))/CElB) 
DEE*DEE+CE1A*DIAG/(1 . 0+ADIFF*ADIFF) 

END  IF 

ELSEIF  (N0N.EQ.4)  THEN 
DEE-DEE+CE2A 
END  IF 
ENDIF 

IF  (K . EQ . 2 . AND . I . GT . 1 )  THEN 

B0N-XE(I-1.J.1)+XE(I.J.2)+XE(I-1,J+1.1)+XECI-1,J,2) 

NON-NINT(RON) 

IF  (NON.EQ.l)  THEN 
DEE*DEE+3*CE2A 
ELSEIF  (NON . EQ . 2)  THEN 
DEE»DEE-2*CE2A 
IF  (XE(I-l.J.l) -GT. .5)  THEN 
DEE-DEE-CE2B 

ADIFF-ABS((XP(I,J)-XP(I-11J+1))/CE1B) 
DEE-DEE+CE1A*DIAG/ ( 1 . 0+ADIFF*ADIFF) 

ELSEIF  (XE(I-1 ,J+1,1).GT. .5)  THEN 
DEE*DEE-CE2B 

ADIFF-ABS((XP(I-1,J)-XP(I,J+1))/CE1B) 
DEE*DEE+CE1A*DIAG/(1 . 0+ADIFF*ADIFF) 

ELSE 

DEE*DEE-2*CE2B 

ADIFF-ABS((XP(I,J)-XP(I-1.J+1))/CE1B) 
DEE«DEE+CE1A*DIAG/(1.0+ADIFF*ADIFF) 
ADIFF-ABS((XP(I-1,J)-XP(I.J+1))/CE1B) 
DEE-DEE+CE1A*DIAG/(1 . 0+ADIFF*ADIFF) 

ENDIF 

ELSEIF  (NON . EQ . 3)  THEN 
DEE-DEE+CE2A 

IF  (XEd-l.J+l.l)  .LT.  .5)  THEN 


DEE00790 

DEE00800 

DEE00810 

DEE00820 

DEE00830 

DEE00840 

DEE00850 

DEE00860 

DEE00870 

DEE00880 

DEE00890 

DEE00900 

DEE00910 

DEE00920 

DEE00930 

DEE00940 

DEE00960 

DEE00960 

DEE00970 

DEE00980 

DEE00990 

DEEOIOOO 

DEEOIOIO 

DEE01020 

DEE01030 

DEE01040 

DEE01050 

DEE01060 

DEE01070 

DEE01080 

DEE01090 

DEEOllOO 

DEEOlllO 

DEE01120 

DEE01130 

DEE01140 

DEE01150 

DEE01160 

DEE01170 

DEE01180 

DEE01190 


DEE-DEE -CE2B 

ADIFF-ABS((XP(I ,J)-XP(I-1 ,  J  +  1))/CE1B) 
DEE-DEE+CE1A*DIAG/ ( 1 . 0+ADIFF*ADIFF) 

ELSEIF  (XE(I-l.J.l).LT. .5)  THEN 
DEE-DEE-CE2B 

ADIFF-ABS ( (XP(I- 1 ,J)-XP(I,J+1))/CE1B) 
DEE-DEE+CE1A*DIAG/(1 . 0+ADIFF*ADIFF) 

ENDIF 

ELSEIF  (N0N.EQ.4)  THEN 
DEE-DEE+CE2A 
ENDIF 
ENDIF 

IF  (K . EQ . 2 . AND . I . LT . NX)  THEN 

RON-XE(I . J, 1) +XE(I+1 . J , 2) +XE(I ,  J+l , 1) +XE( I , J , 2) 
NON-NINT(RON) 

IF  (NON.Eq.l)  THEN 
DEE-DEE+3+CE2A 
ELSEIF  (NON.Eq.2)  THEN 
DEE-DEE- 2+CE2A 
IF  (XE(I.J.l) .GT. .5)  THEN 
DEE-DEE- CE2B 

ADIFF-ABS( (XP(I , J) -XP(I+1,J+1))/CE1B) 
DEE»DEE+CElA*DIAG/( 1 . 0+ADIFF*ADIFF) 
ELSEIF  (XE(I , J+l , 1) . GT. . 6)  THEN 
DEE-DEE-CE2B 

ADIFF-ABS ( (XP( I , J+1)-XP(I+1 , J))/CE1B) 
DEE-DEE+CE1A*DIAG/(1 . O+ADIFF+ADIFF) 

ELSE 

DEE-DEE-2+CE2B 

ADIFF-ABS ( (XP( I ,J)-XP(I+1,J+1))/CE1B) 
DEE-DEE+CE1A*DIAG/(1 . O+ADIFF+ADIFF) 
ADIFF-ABS( (XP(I ,J+1)-XP(I+1,J))/CE1B) 
DEE»DEE+CE1A*DIAG/(1 . O+ADIFF+ADIFF) 

ENDIF 

ELSEIF  (NON.Eq.3)  THEN 
DEE-DEE+CE2A 

IF  (XEU.J+l.l)  .LT.  .5)  THEN 
DEE-DEE-CE2B 

ADIFF-ABS((XP(I. J) -XP (1+1 , J+l) )/CElB) 
DEE-DEE+CE1A*DIAG/(1 . O+ADIFF+ADIFF) 
ELSEIF  (XE(I.J.l) .LT. .5)  THEN 


DEE01200 

DEE01210 

DEE01220 

DEE01230 

DEE01240 

DEE01250 

DEE01260 

DEE01270 

DEE01280 

DEE01290 

DEE01300 

DEE01310 

DEE01320 

DEE01330 

DEE01340 

DEE01350 

DEE01360 

DEE01370 

DEE01380 

DEE01390 

DEE01400 

DEE01410 

DEE01420 

DEE01430 

DEE01440 

DEE01450 

DEE01460 

DEE01470 

DEE01480 

DEE01490 

DEE01500 

DEE01510 

DEE01520 

DEE01530 

DEE01540 

DEE01550 

DEE01560 

DEE01670 

DEE01580 

DEE01590 

DEE01600 


v?  ,»:»  >;y»-  v-w,-.* 


v^-r.y: 


DEE*DEE-CE2B  DEE01610 

ADIFF«ABS( (XP(I . J+l) -XP(I*1 . J) )/CElB)  DEE01620 

DEE«DEE+CE1A*DIAG/(1 . 0+ADIFF*ADIFF)  DEE01630 

END IF  DEE01640 

ELSEIF  (NON . EQ . 4)  THEN  DEE01650 

DEE-DEE+CE2A  DEE01660 

END IF  DEE01670 

ENDIF  DEE01680 

C  CONTRIBUTION  FORM  INHIBITION  OF  PARALLEL  LINES  DEE01690 

IF  (K.EQ.l)  THEN  DEE01700 

DEE*DEE+CE2C* (XE(I-2 ,J,1)+XE(I-1,J,1)+XE(I+1,J , 1) +XE(I+2 , J , 1) )  DEE01710 
ELSE  DEE01720 

DEE*DEE+CE2C* (XE( I , J-2 , 2) +XE(I , J-l , 2) +XE(I ,  J+l , 2) +XE(I , J+2 . 2) )  DEE01730 
ENDIF  DEE01740 

XE( I . J ,K)*HOLD  DEE01750 

END  DEE01760 


4.10  Function  Subprogram  GGUBFS. 

The  function  subprogram  GGUBFS  is  from  the  proprietary  IMSL  Library  and  is  used 
to  generate  pseudorandom  numbers  that  are  independent  and  uniformly  distributed  on 
(0,1).  The  listing  below  should  not  be  reproduced  nor  incorporated  in  any  programs  other 
than  the  present  one  unless  its  use  is  licensed  on  the  system  on  which  such  a  program  is 
being  developed. 


C  IMSL  ROUTINE  NAME  -  GGUBFS  GGU00010 

C  GGU00020 

C . - - - - - - - - - - GGU00030 

C  GGU00040 

C  COMPUTER  -  IBM/SINGLE  GGU00050 

C  GGU00060 

C  LATEST  REVISION  -  JUNE  1.  1980  GGU00070 

C  GGU00080 

C  PURPOSE  -  BASIC  UNIFORM  (0,1)  RANDOM  NUMBER  GENERATOR  -  GGU00090 

C  FUNCTION  FORM  OF  GGUBS  GGU00100 

C  GGU00110 

C  USAGE  -  FUNCTION  GGUBFS  (DSEED)  GGU00120 

C  GGU00130 

C  ARGUMENTS  GGUBFS  -  RESULTANT  DEVIATE.  GGU00140 

C  DSEED  -  INPUT/OUTPUT  DOUBLE  PRECISION  VARIABLE  GGU00150 

C  ASSIGNED  AN  INTEGER  VALUE  IN  THE  GGU00160 

C  EXCLUSIVE  RANGE  (l.DO,  2147483647 , DO) .  GGU00170 

C  DSEED  IS  REPLACED  BY  A  NEW  VALUE  TO  BE  GGU00180 

C  USED  IN  A  SUBSEQUENT  CALL.  GGU00190 

C  GGU00200 

C  PRECISION/HARDWARE  -  SINGLE/ALL  GGU00210 

C  GGU00220 

C  REQD.  IMSL  ROUTINES  -  NONE  REQUIRED  GGU00230 

C  GGU00240 

C  NOTATION  -  INFORMATION  ON  SPECIAL  NOTATION  AND  GGU00250 

C  CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL  GGU00260 

C  INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP  GGU00270 

C  GGU00280 

C  COPYRIGHT  -  1978  BY  IMSL,  INC.  ALL  RIGHTS  RESERVED.  GGU00290 

C  GGU00300 

C  WARRANTY  -  IMSL  WARRANTS  ONLY  THAT  IMSL  TESTING  HAS  BEEN  GGU00310 

C  APPLIED  TO  THIS  CODE.  NO  OTHER  WARRANTY.  GGU00320 

C  EXPRESSED  OR  IMPLIED.  IS  APPLICABLE.  GGU00330 

C  GGU00340 
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c . - . . . . GGU00350 

C  GGU00360 

REAL  FUNCTION  GGUBFS  (DSEED)  GGU00370 

C  SPECIFICATIONS  FOR  ARGUMENTS  GGU00380 

DOUBLE  PRECISION  DSEED  GGU00390 

C  SPECIFICATIONS  FOR  LOCAL  VARIABLES  GGU00400 

DOUBLE  PRECISION  D2P31M.D2P31  GGU00410 

C  D2P3lM=(2**31)  -  1  GGU00420 

C  D2P31  =(2**31) (OR  AN  ADJUSTED  VALUE)  GGU00430 

DATA  D2P31M/2147483647 . DO/  GGU00440 

DATA  D2P31  /2147483648 . DO/  GGU00450 

C  FIRST  EXECUTABLE  STATEMENT  GGU00460 

DSEED  =  DMOD ( 16807 . DO*DSEED , D2P31M)  GGU00470 

GGUBFS  =  DSEED  /  D2P31  GGU00480 

RETURN  GGU00490 

END  GGU00500 
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5  FLIR  EXAMPLES 


The  algorithm  described  in  Sections  2  and  3  and  implemented  by  the  program  of  Section 
4  has  been  applied  to  a  variety  of  FLIR  images  provided  by  the  Advanced  Modeling 
Team  at  NV&EOL.  The  results  of  selected  experiments  are  included  here. 

For  these  experiments,  the  model  parameters  were  set  on  the  basis  of  inspection  of  the 
digitized  FLIR  images  to  determine  attributes  such  as  dynamic  range  and  noise-variance 
and  on  the  basis  of  the  insights  and  interpretations  of  the  model  parameters  described 
in  Section  3. 

In  each  of  the  photographs  in  Appendix  B,  the  upper-left  panel  contains  a  32  x  32 
section  of  the  observed  image.  The  upper-right  panel  contains  the  result  of  fifty  iterations 
of  the  stochastic  relaxation  algorithm,  with  annealing.  The  lower-left  panel  contains  the 
original  observed  image  plus  additional  noise  having  standard  deviation  8.  The  lower- 
right  panel  contains  the  result  of  fifty  iterations  of  the  stochastic  relaxation  algorithm, 
with  annealing,  applied  to  the  noise  corrupted  image. 

The  model  and  program  parameters  for  the  experiments  are  given  in  the  following 
table; _ 


Model 

Program 

Value 

CElA 

10.4 

B 

CElB 

4.0 

•* 

CE2B 

1.66 

0s 

CE2A 

1.0 

£i 

-4 

-3 

-2 

£4 

-1 

£5 

-1 

£e 

0 

PMIN 

40 

PMAX 

238 

MAXDA 

100 

NDA 

100 

For  the  original  observed  images,  the  standard  deviation  SIGMA  of  the  additive  noise 
presumed  to  be  degrading  the  ideal  image  was  set  to  5. 

For  the  images  to  which  noise  was  added,  the  standard  deviation  in  the  restoration 
algorithm  was  set  to  \/25  +  64  =  9.43. 

Eight  figures  are  included  in  Appendix  B. 
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A  COMPLEX  SYSTEMS  WORKING  PAPERS 


During  the  course  of  the  modeling  project,  a  number  of  internal  working  papers  were 
prepared  describing  progress  and  research  plans  for  specific  aspects  of  the  research  effort. 
These  papers  were  not  intended  for  general  distribution.  Nonetheless,  because  of  the 
direct  cooperation  with  the  Advanced  Modeling  Team  at  NV&EOL,  the  working  papers 
were  all  shared  with  the  leaders  of  the  team.  Titles  of  the  working  papers  directly  related 
to  the  image  analysis  problems  at  NV&EOL  include: 

•  An  entropy  approach  to  relaxation  time,  April  1983. 

•  Updating  schemes  for  image  processing,  June  1983. 

•  Parameter  estimation  for  some  Markov  random  fields,  August  1983. 

•  Synthesis  of  partition  patterns,  August  1983. 

•  Synthesis  of  surface  patterns,  August  1983. 

•  A  computer  experiment  with  sweep  areas,  October  1983. 

•  Some  experiments  with  partition,  shape,  and  network  patterns,  October  1983. 


•  Simulating  cold  patterns  is  difficult,  November  1983. 


•  Stochastic  relaxation  for  some  continuous  generator  spaces,  November  1983. 

•  Remarks  on  annealing  schedules,  December  1983. 

•  Recognizing  objects,  March  1984. 

•  Non-localized  generators,  May  1984. 

•  Parameter  estimation  for  Markov  random  fields  with  hidden  variables  and  experi¬ 
ments  with  the  EM  algorithm,  August  1984. 

•  Aspects  of  image  processing,  September  1984. 

•  Software  for  image  processing  experiments,  November  1984. 

•  Preliminaries  to  target  identification  in  IR-pictures,  April  1985. 

•  Recognizing  patterns  in  the  presence  of  nuisance  parameters,  February  1986. 
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